Workshop: Exploring the InsideAirBnB dataset - Prediction

Preamble

# Clear workspace
rm(list=ls()); graphics.off() 

### Load packages
library(tidyverse) # Collection of all the good stuff like dplyr, ggplot2 ect.
library(magrittr) # For extra-piping operators (eg. %<>%)
library(skimr) # For nice data summaries

Load the data

listings <- read_csv('http://data.insideairbnb.com/denmark/hovedstaden/copenhagen/2020-06-26/data/listings.csv.gz')
listings %>% head()

calendar <- read_csv('http://data.insideairbnb.com/denmark/hovedstaden/copenhagen/2020-06-26/data/calendar.csv.gz')
calendar %>% head()

reviews <- read_csv('http://data.insideairbnb.com/denmark/hovedstaden/copenhagen/2020-06-26/data/reviews.csv.gz')
reviews %>% head()

Data munging/tidying

Inspect again…

data <- listings 

data %>% glimpse()

Variable transformations

data %<>%
  mutate(price = price %>% parse_number(),
         cleaning_fee = parse_number(cleaning_fee),
    price_all = price + cleaning_fee) %>%
  filter(number_of_reviews >= 2)

Varriable selection

data %<>% 
  rename(y = price_all) %>%
  select(y, review_scores_rating, neighbourhood_cleansed, accommodates, room_type, bathrooms,is_business_travel_ready,
         number_of_reviews, cancellation_policy, host_is_superhost, host_identity_verified, bedrooms) 

Filtring observations

data %<>% 
  drop_na(y) %>%
  filter(percent_rank(y) <0.95)

Also check the categorical variables for rare types

data %>% count(room_type, sort = TRUE)

data %<>%
  filter(!(room_type %in% c('Shared room', 'Hotel room')))

Misssing data

listings %<>%
    mutate(across(is_character, ~ifelse(.x == "", NA, .x)))

library(VIM)

data %>%
  aggr(numbers = TRUE, prop = c(TRUE, FALSE))

Looks fine, litte missing

EDA

data %>%
  skim()

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             14064     
Number of columns          12        
_______________________              
Column type frequency:               
  character                3         
  logical                  3         
  numeric                  6         
________________________             
Group variables            None      

── Variable type: character ─────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable          n_missing complete_rate   min   max empty n_unique whitespace
1 neighbourhood_cleansed         0             1     5    25     0       11          0
2 room_type                      0             1    12    15     0        2          0
3 cancellation_policy            0             1     8    27     0        4          0

── Variable type: logical ───────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable            n_missing complete_rate  mean count                
1 is_business_travel_ready         0          1    0     FAL: 14064           
2 host_is_superhost                2          1.00 0.143 FAL: 12051, TRU: 2011
3 host_identity_verified           2          1.00 0.448 FAL: 7762, TRU: 6300 

── Variable type: numeric ───────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable        n_missing complete_rate   mean      sd    p0   p25   p50   p75   p100 hist 
1 y                            0         1     974.   365.      100   698   948  1246 1849   ▂▇▇▆▂
2 review_scores_rating        40         0.997  95.5    5.24     20    93    97   100  100   ▁▁▁▁▇
3 accommodates                 0         1       3.17   1.41      1     2     3     4   12   ▇▅▁▁▁
4 bathrooms                    3         1.00    1.05   0.223     0     1     1     1    5.5 ▇▁▁▁▁
5 number_of_reviews            0         1      19.1   30.0       2     5    10    22  637   ▇▁▁▁▁
6 bedrooms                     5         1.00    1.46   0.780     0     1     1     2    8   ▇▃▁▁▁

library(GGally)
data %>% ggpairs() 

Preprocessing

library(tidymodels)

Train and test split

data_split <- initial_split(data, prop = 0.75, strata = y)

data_train <- data_split  %>%  training()
data_test <- data_split %>% testing()

data_recipe <- data_train %>%
  recipe(y ~.) %>%
  step_center(all_numeric(), -all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes()) %>%
  step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE) 

Models & Workflows

Define Models

model_lm <- linear_reg(mode = 'regression') %>%
  set_engine('lm') 

model_xg <- boost_tree(mode = 'regression', 
                       trees = 100,
                       mtry = tune(), 
                       min_n = tune(), 
                       tree_depth = tune(), 
                       learn_rate = tune()
                       ) %>%
  set_engine("xgboost") 

Workflows

workflow_general <- workflow() %>%
  add_recipe(data_recipe) 

workflow_lm <- workflow_general %>%
  add_model(model_lm)

workflow_xg <- workflow_general %>%
  add_model(model_xg)

Hyperparameter Tuning

Resampling

data_resample <- data_train %>% 
  vfold_cv(strata = y,
           v = 3,
           repeats = 2)

Hyperparameter Tuning

tune_xg <-
  tune_grid(
    workflow_xg,
    resamples = data_resample,
    grid = 10
  )

tune_xg %>% autoplot()

best_param_xg <- tune_xg %>% select_best(metric = 'rmse')
best_param_xg

tune_xg %>% show_best(metric = 'rmse', n = 1)

Fit models

Fit models with tuned hyperparameters

workflow_final_xg <- workflow_xg %>%
  finalize_workflow(parameters = best_param_xg)

fit_lm <- workflow_lm %>%
  fit(data_train)

fit_xg <- workflow_final_xg %>%
  fit(data_train)

Predict

pred_collected <- tibble(
  truth = data_test %>% pull(y),
  base = mean(truth),
  lm = fit_lm %>% predict(new_data = data_test) %>% pull(.pred),
  xg = fit_xg %>% predict(new_data = data_test) %>% pull(.pred),
  ) %>% 
  pivot_longer(cols = -truth,
               names_to = 'model',
               values_to = '.pred')

Evaluate

Metrics

pred_collected %>%
  group_by(model) %>%
  rmse(truth = truth, estimate = .pred) %>%
  select(model, .estimate) %>%
  arrange(.estimate)

pred_collected %>%
  group_by(model) %>%
  rsq(truth = truth, estimate = .pred) %>%
  select(model, .estimate) %>%
  arrange(.estimate)

Visuals

pred_collected %>%
  ggplot(aes(x = truth, y = .pred, color = model)) +
  geom_abline(lty = 2, color = "gray80", size = 1.5) +
  geom_point(alpha = 0.5) +
  labs(
    x = "Truth",
    y = "Predicted price",
    color = "Type of model"
  )

Variable Importance

Linear model

fit_lm %>% 
#  pluck(".workflow", 1) %>%   
  pull_workflow_fit() %>% 
  vip::vip(num_features = 10)

XGBoost

fit_xg %>% 
#  pluck(".workflow", 1) %>%   
  pull_workflow_fit() %>% 
  vip::vip(num_features = 10)